eegmp_analysis/tools/[MEG]PLS/MEGPLS_PIPELINE_v2.02b/MEGPLS_COMPONENTS/MEGpipeline_PreprocMEG.m

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Reads raw data with Fieldtrip & applies preprocessing options. %
% Inputs: Raw MEG data compatible with Fieldtrip.                %
% Last modified: Jan. 15, 2014                                   %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Copyright (C) 2013-2014, Michael J. Cheung
%
% This file is a part of the MEG & PLS Pipeline (MEGPLS). For more 
% details, see the documentation included with the software package.
%
% MEGPLS is free software: you can redistribute it and/or modify it under
% the terms of the GNU General Public License version 2 as published by 
% the Free Software Foundation. This program is distributed in the hope 
% that it will be useful, but WITHOUT ANY WARRANTY; without even the 
% implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. 
% See the GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License along
% with this program. If not, you can download the license here: 
% <http://www.gnu.org/licenses/old-licenses/gpl-2.0>.


function MEGpipeline_PreprocMEG(PreprocSettingsMat, PreprocInputMEG)


% Make sure java paths added:
UseProgBar = CheckParforJavaPaths;


% Load .mat files from preprocessingGUI:
LoadSettings = load(PreprocSettingsMat);

time  = LoadSettings.time;
name  = PreprocInputMEG.name;
paths = PreprocInputMEG.paths;

TargetMarkers  = PreprocInputMEG.epoch.TargetMarkers;
IncludedEvents = PreprocInputMEG.epoch.IncludedEvents;
ExcludedEvents = PreprocInputMEG.epoch.ExcludedEvents;
DataFiletype   = PreprocInputMEG.gui.DataFiletype;


% Set up errorlog and diary:
if exist('ErrorLog_PreprocMEG.txt', 'file')
    system('rm ErrorLog_PreprocMEG.txt');
end
if exist('Diary_PreprocMEG.txt', 'file')
    system('rm Diary_PreprocMEG.txt');
end

diary Diary_PreprocMEG.txt
ErrLog = fopen('ErrorLog_PreprocMEG.txt', 'a');



%==============================================%
% EPOCH/READ IN AND PREPROCESS INPUT MEG DATA: %
%==============================================%

NumSubj = length(name.SubjID);

if UseProgBar == 1
    ppm = ParforProgMon('PREPROCESSING MEG DATASETS:  ', NumSubj, 1, 700, 80);
end

parfor s = 1:length(name.SubjID)
    
    % Check if files specified and exist:
    if isempty(paths.DataFullpath{s})
        disp(['Warning: Input data for: ',name.SubjID,' was not specified.']);
        disp('Skipping Subject.')
        continue;
    end
    
    if ~exist(paths.DataFullpath{s}, 'file')
		ErrMsg = sprintf(['ERROR: Input MEG file does not exist:'...
			'\n %s \n\n'], paths.DataFullpath{s});
        
        LogErrorMessage(ErrMsg, 'PreprocMEG');
        disp('ERROR: Input MEG file does not exist. Skipping subject.')
        continue;
    end
    
	[~, ~, FileExt] = fileparts(paths.DataFullpath{s});
	

    
    %--- Read Header and Event info from input datasets: ---%
    %-------------------------------------------------------%
    
    % If reading from raw dataset format:
    if ~strcmp(FileExt, '.mat')
        try
            Hdr = [];
            Hdr = ft_read_header(paths.DataFullpath{s});  % For sampling-rate
        catch
            ErrMsg = sprintf(['ERROR: Dataset header info could not be read into FT:'...
                '\n SKIPPING dataset: %s \n\n'], paths.DataFullpath{s});
            
            LogErrorMessage(ErrMsg, 'PreprocMEG');
            disp('ERROR: Failed to read dataset. Header could not be read into FT.')
            continue;
        end
        
        try
            AllEventInfo = [];
            AllEventInfo = ft_read_event(paths.DataFullpath{s});
        catch
            ErrMsg = sprintf(['ERROR: Dataset event info could not be read by FT:'...
                '\n SKIPPING dataset: %s \n\n'], paths.DataFullpath{s});
            
            LogErrorMessage(ErrMsg, 'PreprocMEG');
            disp('ERROR: Failed to read event info from dataset.')
            continue;
        end
    end
    
    % If reading from FieldTrip .mat file:
    if strcmp(FileExt, '.mat')
        AllEventInfo = [];
        
        DatasetMat = [];
        DatasetMat = LoadFTmat(paths.DataFullpath{s}, 'PreprocMEG');
        if isempty(DatasetMat)
            continue;  % error loading
        end
        
        if isfield(DatasetMat, 'AllEventInfo')
            AllEventInfo = DatasetMat.AllEventInfo;
            
        else
            if ~isempty(IncludedEvents) || ~isempty(ExcludedEvents)
                ErrMsg = sprintf(['ERROR: Inc./Exc. events specified but dataset does not'...
                    'contain event info field. \n As a result, cannot search if specified'...
                    'events are within trials. \n Note: AllEventInfo field is present in FT'...
                    '.mat files created through [MEG]PLS. \n SKIPPING dataset: %s \n\n'], ...
                    paths.DataFullpath{s});
                
                LogErrorMessage(ErrMsg, 'PreprocMEG');
                disp('ERROR: Inc./Exc. events specified, but event info field missing.')
                continue;
            end
        end
    end
    
    
    
    %--- Define trial definition info for datasets: ---%
    %--------------------------------------------------%
    
    % Results in Nx3 trial definition where:
    % (N,1) = BegSample of trial relative to start of raw data.
    % (N,2) = EndSample of trial relative to start of raw data.
    % (N,3) = Offset where negative means trial begins prior to marker, etc.
    
    
    % For raw datasets read in as is (no T=0 events have been specified):
    % Ex: Reading in data already epoched or reading data in as continuous.
    if ~strcmp(FileExt, '.mat') && isempty(TargetMarkers)
        EpochTrialInfo = [];
        trl            = [];
        
        if Hdr.nTrials == 1  % Reading continuous dataset
            trl = zeros(1,3);
            trl(1,1) = 1;
            trl(1,2) = Hdr.nSamples * Hdr.nTrials;
            trl(1,3) = -Hdr.nSamplesPre;
        else                 % Reading dataset already epoched
            trl = zeros(Hdr.nTrials, 3);
            for i = 1:Hdr.nTrials
                trl(i,1) = (i-1)*Hdr.nSamples + 1;
                trl(i,2) = (i  )*Hdr.nSamples    ;
                trl(i,3) = -Hdr.nSamplesPre;
            end
        end
        
        EpochTrialInfo = trl;
    end
    
    
    % For raw datasets being epoched by specified T=0 events:
    % Define trial info for each target T=0 event found.
    if ~strcmp(FileExt, '.mat') && ~isempty(TargetMarkers)
        EpochTrialInfo = [];
        
        cfgDefineTrial         = [];
		cfgDefineTrial         = LoadSettings.FTcfg.DefineTrial;
        cfgDefineTrial.dataset = paths.DataFullpath{s};
        
        for e = 1:size(TargetMarkers, 1)
            cfgDefineTrial.trialdef.eventtype  = TargetMarkers{e,1};
            cfgDefineTrial.trialdef.eventvalue = TargetMarkers{e,2};
            
            try
                CurrentEventInfo  = ft_definetrial(cfgDefineTrial);
                EpochTrialInfo    = [EpochTrialInfo; CurrentEventInfo.trl];
            catch
                ErrMsg = sprintf(['\nWARNING: Could not find this target T=0 event'...
                    'in the current dataset. \n - Dataset: %s \n - Event: %s \n\n'], ...
                    paths.DataFullpath{s}, TargetMarkers{e,3});
                
                LogErrorMessage(ErrMsg, 'PreprocMEG');
            end
        end
        
        % Check for situation that no trials defined (no T=0 events found):
        if isempty(EpochTrialInfo)
            ErrMsg = sprintf(['ERROR: Failed to define trials for current dataset.'...
                '\n None of the target T=0 events were found in the dataset.'...
                '\n %s \n\n'], paths.DataFullpath{s});
            
            LogErrorMessage(ErrMsg, 'PreprocMEG');
            disp('ERROR: Failed to define trials for dataset. Specified T=0 events not found.')
            continue;
        end
        
        % Sort trials in order by samples:
		EpochTrialInfo = sortrows(EpochTrialInfo, 1);  
    end
    
    
    % For FieldTrip .mat datasets:
    if strcmp(FileExt, '.mat')
        EpochTrialInfo = [];
        EpochTrialInfo = DatasetMat.sampleinfo;  % sampleinfo is up-to-date with dataset.
        
        % Get offset info for dataset:
        % Note: The .trl field is not always kept up-to-date by FT like sampleinfo is.
        MissingOffsetInfo = 0;
        
        % If immediate .cfg does not contain .trl or .offset, check if nested in .previous:
        % If both fields are nested, cannot tell which one is more recent:
        if ~isfield(DatasetMat.cfg, 'trl') && ~isfield(DatasetMat.cfg, 'offset')
            GetOffsetCfg = ft_findcfg(DatasetMat.cfg, 'offset');
            GetTrlCfg    = ft_findcfg(DatasetMat.cfg, 'trl');
            
            if ~isempty(GetOffsetCfg) && ~isempty(GetTrlCfg)
                MissingOffsetInfo = 1;
            end
        end
        
        % Check if immediate cfg.offset is present:
        if isfield(DatasetMat.cfg, 'offset')
            if length(DatasetMat.cfg.offset) == 1 || ...
                    isequal(size(EpochTrialInfo, 1), length(DatasetMat.cfg.offset))
                EpochTrialInfo(:,3) = DatasetMat.cfg.offset;
            else
                MissingOffsetInfo = 1;
            end
          
        % Otherwise, check for .trl field:
        else
            GetTrlCfg = ft_findcfg(DatasetMat.cfg, 'trl');
            
            if isequal(EpochTrialInfo(:,[1:2]), GetTrlCfg(:,[1:2]))  % If .trls up-to-date
                EpochTrialInfo(:,3) = GetTrlCfg(:,3);
                
            else  % If .trls not up-to-date:
                for t = 1:size(EpochTrialInfo, 1)  % Get offset from respective matching trial
                    MatchTrial = find(ismember(GetTrlCfg(:,[1:2]), EpochTrialInfo(t,[1:2]), 'rows'));
                    
                    if ~isempty(MatchTrial) && length(MatchTrial) == 1
                        EpochTrialInfo(t,3) = GetTrlCfg(MatchTrial, 3);
                    else
                        MissingOffsetInfo = 1;
                    end
                end
            end
        end
        
        % Display warnings:
        if MissingOffsetInfo == 1
            if ~isempty(IncludedEvents) || ~isempty(ExcludedEvents)
                ErrMsg = sprintf(['ERROR: Cannot determine trial offset info for dataset.'...
                    '\n As a result, cannot search if inc./exc. events are within trials.'...
                    '\n To fix, use "ft_redefinetrial" on dataset to define cfg.trl field.'...
                    '\n SKIPPING dataset: %s \n\n'], paths.DataFullpath{s});
                
                LogErrorMessage(ErrMsg, 'PreprocMEG');
                disp('ERROR: Cannot search for inc./exc. events for dataset. See ErrLog.')
                continue;
            end
            
            ErrMsg = sprintf(['WARNING: Cannot determine trial offset info for dataset.'...
                '\n As a result, could not check if trials exceeded specified overlap.'...
                '\n All trials will be kept included regardless of overlap amount.'...
                '\n To fix, use "ft_redefinetrial" on dataset to define cfg.trl field.'...
                '\n Dataset: %s \n\n'], paths.DataFullpath{s});
            
            LogErrorMessage(ErrMsg, 'PreprocMEG');
        end
    end  % If .mat
    
    
        
    %--- Remove trials that exceed overlap threshold: ---%
    %----------------------------------------------------%
    
    % Note: The "hdr" field in FT .mat files is not kept up-to-date.
    % - It represents the header of the ORIGINAL .ds file read into FT.
    % - Therefore, if reading in FT .mat, get sampling rate from fsample.
    
    if ~strcmp(FileExt, '.mat')
        SamplingRate = Hdr.Fs;  % Get from .ds
    elseif strcmp(FileExt, '.mat')
        SamplingRate = DatasetMat.fsample;  % Get from .mat structure
    end
    ThresholdSamples = round(time.OverlapThresh * SamplingRate);
    
    RejectTrials = [];
    for n = 2:size(EpochTrialInfo, 1)
        PrevTrialEnd   = EpochTrialInfo(n-1, 2);
        CurrTrialStart = EpochTrialInfo(n,1);
        SampleOverlap  = PrevTrialEnd - CurrTrialStart;
        
        if SampleOverlap < 0
            continue;  % No overlap with previous trial.
        else
            if SampleOverlap > ThresholdSamples
                RejectTrials = [RejectTrials, n-1, n];  % Reject both overlap trials.
            end
        end
    end
    
    RejectTrials = unique(RejectTrials);
    EpochTrialInfo(RejectTrials,:) = [];
    
    
    
    %--- Isolate for trials with included events in search-range: ---%
    %----------------------------------------------------------------%
    
    if ~isempty(IncludedEvents)
        
        % Get sample info of when included events occurred:
        IncEventSamples = [];
        for e = 1:size(IncludedEvents, 1);
            EventType = IncludedEvents{e,1};
            EventInfo = AllEventInfo(ismember({AllEventInfo.type}, EventType));
            
            EventValue = IncludedEvents{e,2};
            if ~isempty(EventValue)
                EventInfo = EventInfo(ismember([EventInfo.value], EventValue));
            end
            
            IncEventSamples = [IncEventSamples, EventInfo.sample];
            IncEventSamples = sort(IncEventSamples);
        end
        
        % Get time-range of data (in samples) to search for IncludedEvents:
        SampleSearchRange = [];
        for n = 1:size(EpochTrialInfo, 1)
            TrialBegSample = EpochTrialInfo(n,1);
            TrialEndSample = EpochTrialInfo(n,2);
            TrialOffset    = EpochTrialInfo(n,3);
            
            StartSearchSample = TrialBegSample - TrialOffset + ...
                (round(time.IncEventStart * SamplingRate));
            
            EndSearchSample   = TrialBegSample - TrialOffset + ...
                (round(time.IncEventEnd * SamplingRate));
            
            if StartSearchSample < TrialBegSample
                StartSearchSample = TrialBegSample;
            end
            if EndSearchSample > TrialEndSample
                EndSearchSample = TrialEndSample;
            end
            
            SampleSearchRange(n,:) = [StartSearchSample:EndSearchSample];
        end
        
        % Only keep trials where included events found in search-range:
        RejectTrials = [];
        for n = 1:size(EpochTrialInfo, 1)
            if max(ismember(SampleSearchRange(n,:), IncEventSamples)) ~= 1
                RejectTrials = [RejectTrials, n];
            end
        end
        
        RejectTrials = unique(RejectTrials);
        EpochTrialInfo(RejectTrials,:) = [];
    end
    
    
    
    %--- Exclude trials with excluded events in search-range: ---%
    %------------------------------------------------------------%
    
    if ~isempty(ExcludedEvents)
        
        % Get sample info of when excluded events occurred:
        ExcEventSamples = [];
        for e = 1:size(ExcludedEvents, 1);
            EventType = ExcludedEvents{e,1};
            EventInfo = AllEventInfo(ismember({AllEventInfo.type}, EventType));
            
            EventValue = ExcludedEvents{e,2};
            if ~isempty(EventValue)
                EventInfo = EventInfo(ismember([EventInfo.value], EventValue));
            end
            
            ExcEventSamples = [ExcEventSamples, EventInfo.sample];
            ExcEventSamples = sort(ExcEventSamples);
        end
        
        % Get time-range of data (in samples) to search for ExcludedEvents:
        SampleSearchRange = [];
        for n = 1:size(EpochTrialInfo, 1)
            TrialBegSample = EpochTrialInfo(n,1);
            TrialEndSample = EpochTrialInfo(n,2);
            TrialOffset    = EpochTrialInfo(n,3);
            
            StartSearchSample = TrialBegSample - TrialOffset + ...
                (round(time.ExcEventStart * SamplingRate));
            
            EndSearchSample   = TrialBegSample - TrialOffset + ...
                (round(time.ExcEventEnd * SamplingRate));
            
            if StartSearchSample < TrialBegSample
                StartSearchSample = TrialBegSample;
            end
            if EndSearchSample > TrialEndSample
                EndSearchSample = TrialEndSample;
            end
            
            SampleSearchRange(n,:) = [StartSearchSample:EndSearchSample];
        end
        
        % Reject trials where excluded events found in search-range:
        RejectTrials = [];
        for n = 1:size(EpochTrialInfo, 1)
            if max(ismember(SampleSearchRange(n,:), ExcEventSamples)) == 1
                RejectTrials = [RejectTrials, n];
            end
        end
        
        RejectTrials = unique(RejectTrials);
        EpochTrialInfo(RejectTrials,:) = [];
    end
    
    
    %--- Epoch/Read in data and apply preprocessing: ---%
    %---------------------------------------------------%
    
    if isempty(EpochTrialInfo)
        ErrMsg = sprintf(['ERROR: Current event settings result in no trials for dataset:'...
            '\n %s \n\n'], paths.DataFullpath{s});
        
        LogErrorMessage(ErrMsg, 'PreprocMEG');
        disp('ERROR: Event parameters result in no trials for current dataset.')
        continue;
    end
    

    % If raw dataset format:
    if ~strcmp(FileExt, '.mat')
        cfgPreproc         = [];
        cfgPreproc         = LoadSettings.FTcfg.Preproc;
        cfgPreproc.dataset = paths.DataFullpath{s};
        cfgPreproc.channel = 'all';
        cfgPreproc.method  = 'trial';
        cfgPreproc.trl     = EpochTrialInfo;
        
        if Hdr.nTrials == 1  % If raw dataset is continuous
            cfgPreproc.padtype = 'data';
        else                 % If raw dataset already epoched
            cfgPreproc.padtype = 'mirror';
        end
        
        MEGdata = [];
        MEGdata = ft_preprocessing(cfgPreproc);
    
        
    % If FieldTrip .mat:
    elseif strcmp(FileExt, '.mat')
        cfgPreproc         = [];
        cfgPreproc         = LoadSettings.FTcfg.Preproc;
        cfgPreproc.channel = 'all';
        cfgPreproc.method  = 'trial';
        cfgPreproc.padtype = 'mirror';
        
        DatasetMat           = [];  % Free memory
        cfgPreproc.inputfile = paths.DataFullpath{s};
        
        MEGdata = [];
        MEGdata = ft_preprocessing(cfgPreproc);
        
        if MissingOffsetInfo == 0
            cfgRedefine     = [];
            cfgRedefine.trl = EpochTrialInfo;
            MEGdata = ft_redefinetrial(cfgRedefine, MEGdata);
        end
    end
    
    
    % Run CTF gradient-correction if specified:
    cfgGradCorrectCTF  = [];
    cfgGradCorrectCTF  = LoadSettings.FTcfg.GradCorrectCTF;
    
    if ~strcmp(cfgGradCorrectCTF, 'none')
        MEGdata = ft_denoise_synthetic(cfgGradCorrectCTF, MEGdata);
    end
    
    
    % Add event info into MEGdata structure in case of future preprocessing:
    if ~isfield(MEGdata, 'AllEventInfo')
        MEGdata.AllEventInfo = AllEventInfo;
    end
    
    
    
    %--- Save preprocessed data: ---%
    %-------------------------------%
    
    GroupFolder = ['GroupID_',name.CurrentGroupID];
    CondID      = name.CurrentCondID;
    
    OutputFullpath = ...
        [paths.DataID,'/',GroupFolder,'/',CondID,'/',name.SubjID{s},'_PreprocMEG.mat'];
    
    CheckSavePath(OutputFullpath, 'PreprocMEG');
    ParforSaveMEGdata(OutputFullpath, MEGdata);
    
    MEGdata        = [];  % Free memory
    EpochTrialInfo = [];
    AllEventInfo   = [];
    Hdr            = [];
    
    if UseProgBar == 1 && mod(s, 1) == 0
        ppm.increment();  % move up progress bar
    end
    
end  % Subj
if UseProgBar == 1
    ppm.delete();
end



%=========================%
% CHECK FOR OUTPUT FILES: %
%=========================%

for s = 1:length(name.SubjID)
    GroupFolder = ['GroupID_',name.CurrentGroupID];
    CondID      = name.CurrentCondID;
    
    OutputFullpath = ...
        [paths.DataID,'/',GroupFolder,'/',CondID,'/',name.SubjID{s},'_PreprocMEG.mat'];
    
    if ~exist(OutputFullpath, 'file')
        fprintf(ErrLog, ['ERROR: Missing MEG preprocessed output file for:'...
            '\n %s \n\n'], OutputFullpath);
    end
end



%=================%

if exist([pwd,'/ErrorLog_PreprocMEG.txt'], 'file')
    LogCheck = dir('ErrorLog_PreprocMEG.txt');
    if LogCheck.bytes ~= 0  % File not empty
        open('ErrorLog_PreprocMEG.txt');
    else
        delete('ErrorLog_PreprocMEG.txt');
    end
end

fclose(ErrLog);
diary off

end



function ParforSaveMEGdata(OutputPath, MEGdata)
    MEGdata;
    save(OutputPath, 'MEGdata');
end